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^ . ABSTRACT 



We solve the nonhnear evolution of pressureless, irrotational density 
fluctuations in a perturbed Robert son- Walker spacetime using a new Lagrangian 
method based on the velocity gradient and gravity gradient tensors. Borrowing 
^ I results from general relativity, we obtain a set of Newtonian ordinary differential 

I equations for these quantities following a given mass element. Using these 

Lagrangian fluid equations we prove the following collapse theorem: A mass 
I element whose density exceeds the cosmic mean at high redshift collapses to 

^ ' infinite density at least as fast as a uniform spherical perturbation with the 

same initial density and velocity divergence. Velocity shear invariably speeds 
collapse — the spherical tophat perturbation, having zero shear, is the slowest 
configuration to collapse for a given initial density and growth rate. Two 



2 ' corollaries follow: (1) Initial density maxima are not generally the sites where 



O ■ 

^ ■ collapse first occurs. The initial velocity shear (or tidal gravity field) also is 

important in determining the collapse time. (2) Initially underdense regions 
^ ■ undergo collapse if the shear is sufficiently large. 

' If the magnetic part of the Weyl tensor vanishes, the nonlinear evolution is 

described purely locally by these equations. This condition is exact for highly 
symmetrical perturbations (e.g., with planar, cylindrical, or spherical symmetry) 
and may be a good approximation in many other circumstances. Assuming 
the vanishing of the magnetic part of the Weyl tensor we compute the exact 
nonlinear gravitational evolution of cold matter. We find that 56% of initially 
underdense regions collapse in an Einstein-de Sitter universe for a homogeneous 
and isotropic random field. We also show that, given this assumption, the final 
stage of collapse is generically two-dimensional, leading to strongly prolate 
filaments rather than Zel'dovich pancakes. While this result may explain the 
prevalence of filamentary collapses in some N-body simulations, it is not true 
in general, suggesting that the magnetic part of the Weyl tensor does not 
necessarily vanish in the Newtonian limit. 
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Subject headings: cosmology: theory — large-scale structure of the universe — 
gravitation 



1. Introduction 



According to the gravitational instability theory, cosmic structure developed from the 
amplification of small density fluctuations generated in the early universe. During the 
linear stage the amplitude is small and the evolution is described simply by the linearized 
equations of motion (gravitational field equations plus Boltzmann or fluid equations for the 
matter and radiation components). If pressure gradients are not too large, self-gravitating 
fluctuations eventually become nonlinear. The difficulty of following the nonlinear evolution 
analytically has led cosmologists to resort to gravitational N-body simulations for studying 
the formation of galaxies and large-scale structure (e.g., [Bertschinger 199"l| and references 
therein). However, these simulations should be accompanied by analytical models enabling 
one to check and comprehend the results. 

There are two nonlinear models that have been widely used to describe nonlinear 
gravitational evolution in cosmology. The first is the spherical model, for which a simple 
exact solution exists ( [Tolman 1934 ; Pondi 1947 ; Partridge fc Peebles 1967 ; Gunn fc Gott 



|1972| ; Peebles 19801 , §19). For pressureless spherical collapse one follows the trajectories of 
spherical mass shells having fixed enclosed mass. This model is generally applied to the 
formation of dense objects, the logic being that the first structures to collapse are high 
density peaks, with a nearly spherical distribution of matter surrounding them (pardeen et 
|al. 19861 ). 

The second model is the kinematical one of Zel'dovich (1970). The key idea here is to 
extrapolate the linear evolution of peculiar velocity beyond the linear regime (for a review 
see [Shandarin fc Zel'dovich 1989|) . The Zel'dovich method is also Lagrangian but is not 



restricted to spherical symmetry. Indeed, for plane-parallel perturbations (but not for two- 
or three-dimensional ones) it is exact for particles whose trajectories have not intersected 
others. 

The Zel'dovich approximation (as it is known in multi-dimensions) has provided a 
great deal of insight into the initial nonlinear evolution of density fluctuations. Recently 
it has inspired a number of developments concerning the evolution of the velocity field 
for irrotational flows. Nusser et al. (1991) showed how to compute the density from the 
Eulerian velocity field in the Zel'dovich approximation; their work was extended to second 
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order in Lagrangian perturbation theory by Gramann (1993). Lagrangian perturbation 
theory has also been developed and applied by Moutarde et al. (1991) and Bouchet et al. 
(1992). A broader perturbative framework has been developed by Giavalisco et al. (1993). 
The Zel'dovich approximation has also been used to estimate the probability distribution of 
density in the nonlinear regime ([Kofman 1991| ; padmanabhan fc Subramanian 1993| ; [Kofman 



|et al. 1994|) . Bernardeau (1992) has used a somewhat different statistical summation 



method to relate the density and velocity divergence in the nonlinear regime. 

The limitations of the spherical model and Zel'dovich approximation are well known. 
The first applies only with spherical symmetry while the second gives inaccurate trajectories 
beyond linear perturbation theory. Various authors have therefore suggested modifications 
of the Lagrangian approach to increase its accuracy. Buchert (1992) has investigated 
the linear and nonlinear limits of the Zel'dovich approximation. He and Barrow & Saich 
(1993) have also studied the effects of vorticity. Gurbatov, Saichev, & Shandarin (1980) 
and Kofman, Pogosyan, & Shandarin (1990) have proposed adding viscosity to prevent 
the intersection of fluid trajectories. Lachieze-Rey (1993a,b) has recently developed 
elegant matrix methods for extending the Zel'dovich approximation beyond linear order by 
following the evolution of the deformation tensor, the Jacobian of the transformation from 
Lagrangian to Eulerian coordinates. Matarrese et al. (1992) have suggested the interesting 
approximation of freezing the initial fluid streamlines using the linear velocity field, in effect 
treating the velocity potential as constant. Brainerd, Scherrer, & Villumsen (1993) and 
Bagla & Padmanabhan (1994) have suggested instead treating the gravitational potential 
as constant. While these approaches have provided valuable insights about nonlinear 
evolution, they are all based on approximations to the evolution of the gravitational field. 

General relativity can contribute to a solution by providing evolution equations for 
the gravitational field as opposed to the nonlocal static Poisson equation underlying the 
Newtonian approach. It has long been known among relativists that the Einstein equations 
and the Bianchi identities provide both dynamical equations and initial-value constraints 
for the gravitational field, whether specified by the metric or by the Ricci and Weyl tensors 



(see [Ellis 1971| for an excellent overview). Moreover, the relativistic equations naturally 
lead to a closed set of equations for the fluid variables and gravitational fields specified 
by the Weyl tensor. The Weyl tensor may be decomposed into electric and magnetic 
parts; in the Newtonian limit the electric part represents the gravitational tidal field. The 
Bianchi identities provide an exact Lagrangian evolution equation for the Weyl tensor and, 
therefore, for the tidal field in the Newtonian limit. Remarkably, when the magnetic part 
of the Weyl tensor vanishes the evolution of the tidal field is purely local, allowing it to be 
integrated together with the density and the velocity gradient tensor independently for each 
mass element. 



-4- 



The fully relativistic Lagrangian equations of motion have been used to describe the 
linear evolution of cosmological perturbations ( [Hawking 1966| ; [Ellis fc Bruni 1989| ; pwang 



fc Vishniac 199(J|; Bruni, Dunsby, & Ellis 1992). Their application to nonlinear evolution 



of cosmological fluctuations was first made by Barnes & Rowlingson (1989). A major step 
applying these methods to large-scale structure was made by Matarrese, Pantano, & Saez 
(1993), who showed that exact solutions exist prior to orbit-crossing for spherical and 
plane-parallel perturbations, and that these solutions are identical to the afore-mentioned 
spherical model and Zel'dovich solution. They also demonstrated a local mapping between 
Lagrangian and Eulerian coordinates. The Lagrangian fluid method in general relativity 
has been developed further by Croudace et al. (1993), who introduced an approximation 
method based on the Hamilton- Jacobi equation. 

The current paper presents the Lagrangian fluid approach in its Newtonian framework 
and applies it to the general problem of nonlinear gravitational collapse beginning from 
small perturbations of a Friedmann-Robertson- Walker spacetime. In §2 we present the 
Lagrangian fluid equations and their solution in second-order perturbation theory. Although 
these equations have been given by Matarrese et al. (1993), our approach emphasizes the 
Newtonian interpretation in contrast with their general relativistic one. In §3 we present a 
collapse theorem giving a lower bound on the collapse time of an overdense perturbation. 
The full equations are integrated assuming vanishing magnetic part of the Weyl tensor 
in §4, where we present results for both overdense and underdense initial perturbations. 
Conclusions are given in §5. 



2. Nonlinear Evolution of Density, Velocity Gradients, and Tides 



The evolution of a pressureless fluid may be described either by the Euler equations, by 
geodesic equations for individual mass elements, or by Lagrangian fluid equations. Although 
our goal is a closed system of Lagrangian fluid equations, we begin with the nonrelativistic 
Eulerian approach because of its greater familiarity. In a perturbed Robertson- Walker 
universe with expansion scale factor a(r) the Euler equations take the following form 
( [Bertschinger 19"^^ : 

^ + V-[{l + 6)v] = 0, (la) 

— + (v-\/)v = — {/-V0, lb 
ar ^ ^ a 

= AnGa^p = ^ QoHoa'^S . (Ic) 
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The mass density is p = p + 5p = p(t)(1 + 5) and the pecuhar velocity is {7 = dx/dr where 
r is conformal time {dr = dt/a, with dots denoting derivatives with respect to r) and x is 
the comoving position with corresponding gradient V = e*Vj = e^d/dx^ with orthonormal 
basis vectors: Cj ■ Cj = 5ij. Throughout this paper we assume the background space to be 
flat on scales of interest. Equations (1) are valid for a nonrelativistic pressureless fluid on 
scales much smaller than the Hubble distance cH~^. In adopting a Newtonian approach we 
neglect gravitational waves and gravomagnetism (the gravitational effect of vorticity). 

We will rewrite the fluid equations in terms of the Lagrangian time derivative following 
a fixed fluid element, d/dr = d/dr + v ■ V. Each Lagrangian fluid element has a trajectory 
x{q,T) where g is a Lagrangian coordinate so that d/dr is the time derivative at fixed q 
while d/dr is taken at fixed x. Converting the time derivatives from Eulerian to Lagrangian, 
equation ( p!a| ) becomes 

^ + (1 + 5)^ = 0, 9 = W-v. (2) 
dr 

From equation (0) we see that we require an evolution equation for the velocity 
divergence (or expansion scalar). For this purpose we consider the velocity gradient (or 
rate-of-strain) tensor. We decompose this tensor into its trace and traceless symmetric and 
antisymmetric parts in the usual way: 

Vjfj = - 9 6ij + (Ty + tOij , cTjj = Cji , Uij = eijku'^ = —iLjji , (3) 

where 2 tl; = V x -y is the vorticity and eijk is the three-dimensional Levi-Civita symbol. 
The unperturbed Hubble flow is absent from 6 because v is the peculiar velocity. Taking 
the divergence of equation ([T^), using equations (|T^ and (^), and converting the time 
derivative from Eulerian to Lagrangian, one obtains the Raychaudhuri equation: 

^ + ^e + -e'' + - 2^2 = -^T^GpaH , (4) 

dr a 3 

where a;^ = u'^uji. 

Next we require evolution equations for the shear cjjj and vorticity. These follow from 
the traceless symmetric and antisymmetric parts of the gradient of equation ( |lb|) . The curl 
gives 

'^ + ^u^ + leuj^-a%uj^ = Q. (5) 
dr a 6 

The traceless symmetric part of the gradient of equation dTH) gives 



daij d 2 
'd^ 
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where Eij = VjVj0 — (1/3) Sij is the gravitational tidal field, known in general relativity 
as the electric part of the Weyl tensor in the fluid frame ( |Hawking 1966| ; [Ellis 197T| ). 



Equations (H) and (11)-® give a nearly closed system of equations for the mass density 
and the velocity gradient (represented by the expansion, shear, and vorticity) following a 
Lagrangian fluid element. However, the gravitational tidal fleld remains. In the Newtonian 
framework, Eij is found after solving the Poisson equation for (p. Because this is a nonlocal 
operation involving the mass distribution everywhere, it would appear to be impossible to 
obtain a set of Lagrangian equations for one mass element. 

This conclusion is false. In general relativity, the fleld equations are equivalent to a 
set of evolution equations plus constraints. The Poisson equation corresponds to one of 
the constraints, the ADM energy constraint equation. However, there exist evolution and 
constraint equations for all of the gravitational flelds, including the tidal tensor Eij. It is 
possible to obtain an evolution equation for Eij. However, this equation is not necessarily 
local. 

Ellis (1971) gives a clear discussion of the Lagrangian evolution equations and 
constraints for gravitational flelds in general relativity. Equations of motion for the Weyl 
tensor (giving that part of the gravitational fleld not produced by local sources) follow from 
the Bianchi identities rather than the Einstein equations. Assuming that the matter has 
nonrelativistic peculiar velocities and pressure, the evolution equation for the electric part 
of the Weyl tensor ( [Elhs 197T| ; [Matarrese et al. 1993| ; Proudace et al. 1993|) , expressed in 



comoving coordinates, is 

+ ^ + BE,, + 6,, a'^Eki - 3a\,E,)k + e'' uEj)kOOi - V, e'^' ,,Hjy = -AnGpa^a,, . (7) 



Parentheses around a pair of subscripts indicates symmetrization: a^^-Ej-^^ = 

I {a'^iEjk + jEik). The fully antisymmetric Levi-Civita tensor is eijk (with ei23 = +1 in 

Cartesian coordinates). The magnetic part of the Weyl tensor is denoted Hij. We have not 

assumed that density fluctuations or the velocity gradient are small; the total mass density 

p (as opposed to the fluctuation 5p) appears on the right-hand side. As in the preceding 

equations, the time derivative is Lagrangian, so that each fluid element q is assigned a value 

Eij{q,T). 

Aside from the term involving Hij, equation (|^) is purely local. Thus, if Hij = 0, we 
have obtained a closed set of local Lagrangian equations for the nonlinear evolution of 
pressureless irrotational matter. This fact was flrst noted by Barnes & Rowlingson (1989) 
and applied in cosmology by Matarrese et al. (1993, 1994). If H^j = and if equation (|l^) 
is applied as a constraint at some initial time, then integrating equation (^ for Eij gives 
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exactly the same results (prior to orbit-crossing) as if Eij were computed from the solution 
of the Poisson equation at the final time. 

Do local evolution equations for the gravitational tidal field violate causality? No, 
because they follow directly from general relativity, which is manifestly causal. The 
resolution of this paradox is the fact that on account of the Bianchi identities the Einstein 
field equations contain information about the evolution of the matter: the contracted 
Bianchi identities imply local energy-momentum conservation. Causality is preserved as 
long as the motion of the matter is such as to generate no gravitational radiation. 

We shall not write down the evolution equation for Hij here; it is given by Ellis (1971). 
It is clear that if the gradient of H^j (which appears in eq. 0) is nonzero, local evolution 
breaks down and the evolution of the tide and velocity variables becomes much more 
complicated. In this paper we shall either restrict ourselves to proving results that are valid 
independently of the value of Hij, or we shall assume that it can be neglected. While the 
latter choice is not entirely satisfactory, it may be regarded as an approximation whose 
validity is to be established later. By making this approximation we can easily integrate 
the evolution equations and compare the results against N-body or other fully numerical 
solutions. Thus, we can make strong predictions about gravitational collapse that can be 
used to assess the importance of the magnetic part of the Weyl tensor even if we lack a 
good Newtonian interpretation. 

Equation (^) shows that if = initially, it remains zero. This statement remains 
valid for individual coUisionless mass elements even after trajectories intersect. In the 
following we assume that primeval vorticity is negligible so that a; = 0. 

We will see below that for growing-mode fluctuations in the linear regime with = 0, 
Eij oc (Jij and therefore the eigenvectors of these matrices are aligned. If we choose our 
coordinate axes to be aligned with these eigenvectors initially, then equations @ and (^) 
are diagonal, implying that the eigenvectors remain fixed in direction ( Parnes fc Rowlingson 



1989|) . Thus, each of the tensor equations (y) and (0) may be reduced to equations of 
motion for two of the eigenvalues (with the sum of eigenvalues vanishing for a traceless 
matrix). For convenience we write the tensors in the following form: 



<^ij = ^ ^ <5y (a) 5 = ^ Gpa'^ e (1 + 5) Qij{P) , (8) 



where we have introduced new scalars c < 0, e > 0, a, and (3, and a one-parameter traceless 
quadrupole matrix 



[Qijia)] = diag 



a + 27r\ /a — 27r\ /a 



cosl^J,cos^^j,cos^3 



(9) 
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In the following we will refer to cr and e as the shear and tide scalars, respectively. Similarly, 
a and (3 are called the shear and tide angles, respectively. These angles give the ratios of 
eigenvalues (or axis ratios) of the shear and tidal tensors. 

This matrix representation is convenient because all possible eigenvalues are obtained 
by qQij{a) with q G [0, oo) (or —q in the same range) and a G [0,7r]. Although Qij{a) 
is periodic with period 6tt rather than 2tt, Qij{a ± 2tt) and Qij{—a) differ from Qij{a) 
only by permutation of the coordinates. The reduction of the angular range by a factor 
of six corresponds to ordering the eigenvalues so that Qu < Q22 < Qss- Two of the 
eigenvalues are negative for a G [0,7r/2) while two are positive for a G (7r/2,-7r]. The sum 
of the eigenvalues vanishes. The square and determinant of Qij are simple: Q^^Qij = |, 
det [Qij (a)] = \ cos a. We will also need the following differential and product identities: 

dQijia) = ^Qij (^a + ^) da , 2Qik{a)QkjiP) = cos f— j ^ij + Qii(-« - P) ■ (10) 



Substituting equations (||) into equations (|) and (|^) and using equations 



and 



(p!0|), we get the following equations of motion for the shear and tide scalars and angles: 



da a 1 , „ 
— + -a + -ai2e 
ar a 6 



a cos a) = —ATxGpa? e (1 + 5) cos 



a 



13' 



da 
dr 



de 

dr 

d(3 
dr 



asm a 



ae cos 



+ 3a sin 




(11) 
(12) 
(13) 
(14) 



Equations (p!3| ) and (|T^) have assumed Hij = 0. Equations (|Tl])-(|T^) have singular solutions 
with 6, —6, — cr, and e diverging for finite r. When approaching a singularity by numerical 
integration we change the independent variable from t to x = ln(l + 6) with dx = —Odr 
and then integrate r(x) and {6/6){x), {a/9){x), and {e/6){x) in place of 5(r), 0{r), cr(r), 
and e(r). 

With zero vorticity and pressure, equations (H), (||) (with a^^aij becoming |cr^), and 
(p!T|)-([T^) completely determine the evolution of the density, expansion, shear, and tides 
until fluid trajectories intersect. We assume that at early times the perturbations were 
small, 5^ <^ 1, and i7 ^ 1 (this is generally valid even if Q ^ 1 today). The solution of 
the linearized equations then gives S = (5oa(r) and 6 = —6. However, we must also specify 
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the initial eigenvalues of the shear and tidal tensors. For gravitationally-induced linear 
perturbations these two tensors are proportional to each other, as one may see by linearizing 
equation (6). Thus, altogether we need three constants to specify a growing-mode solution. 
We take these to be the the linear density fluctuation 60, tide (or shear) scalar eo, and 
tide (or shear) angle ao- These quantities are simply related to the double gradient of the 
gravitational potential: 



Qq Hq 



1 



SoSij + eo Qij{ao) 



(15) 



where we have assumed that coordinate axes have been oriented locally to coincide with 
the eigenvectors of VjVj0. 



In terms of the perturbation constants {Sq, eo, ao) the second-order solution for small a 



IS 



6 = a6o + ^a\6l + 6,) , 9 = ~a {5, + aO^) , 0^ = ^ (l3 5l + 8e; 



a = —a [eQ + aai) 

e = a(eo + aei) , ei 



cTi = — ( 26eo5o - 5eo cos ao 



13 
21 



a = ao + - aeo sin ao 
13 



(eoSo — el cos ao) , P = ao + — aeo sin ao 



(16a) 
(16b) 
(16c) 



In the linear regime the expansion scalar is proportional to minus the density 
perturbation and the shear tensor is proportional to minus the tidal tensor, implying 
ViVj oc — ViVj0. This proportionahty occurs because for growing-mode perturbations, the 
peculiar velocity is proportional to the gravity vector — V0. Thus, the initial conditions 
may be specified by the eigenvalues of either ViVj0 or — Vjfj, as they are equivalent. 
However, the proportionality of the shear and tide is broken in second-order perturbation 
theory. 



3. Collapse Theorem 



We note from equations (H) and (Q) that shear (nonzero a) increases the rate of 
growth of the convergence —6, thereby increasing the rate of growth of density fluctuations. 
Vorticity (nonzero u) inhibits this process. With zero shear and vorticity the velocity 
gradient tensor is isotropic (eq. 0), corresponding to uniform spherical collapse with radial 
motions toward some center. This suggests that uniform spherical perturbations grow more 
slowly than others. 
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To extend this argument into a theorem we combine equations and into a 
second-order differential equation for S{t): 

S+^S = ^-^^ + {l + S)(la'-2u;' + 4nGpa'6) . (17) 
a 6 1 + \6 J 

This equation is exact for pressureless matter, regardless of whether the magnetic part of 
the Weyl tensor vanishes. We see that if a; = 0, the growth rate of 5 is minimized for a = 0. 
This conclusion holds independently of assumptions about the evolution of the shear or 
tides and is due to the simple geometrical fact that shear increases the rate of growth of the 
convergence of fluid streamlines. 

When (T = = 0, equation ([T7|) reduces to the exact equation for the evolution of 
the mean density in the spherical model. In other words, this equation is equivalent to 
the equation of motion for a spherical shell of proper radius r(t) enclosing a fixed mass M 
related to 5 by 1 + 5 = 3M/(47rpr'^). From the exact solution it is known that spherically 
symmetric growing density fluctuations collapse if 5j = Oj^o > (3/5)(i7j~^ — 1) where 5j 
and f2j refer to the density fluctuation and cosmic density parameter at initial expansion 
parameter ([Peebles 1980| , §19E). If = i7j = 1, the collapse of a spherical overdense 
perturbation occurs at a = air^ = ac obeying the well-known relation 

5 / 2 \ 

a;i(5o>0,eo = 0) = - — 6o = C,6o, (18) 



3 VStt 



with Ci = 0.592954. . .. In other words, if linear theory would imply 5 = 1 at a = 1, then 
the density of a uniform spherical perturbation becomes infinite at a = Cf^ = 1.68647. . .. 

Based on the above considerations, we are led to state the following Collapse Theorem 
for overdense growing fluctuations: A pressureless, irrotational mass element with initial 
density fluctuation 6i > {3/5){Ql^ — 1) collapses at least as fast as a uniform spherical 
perturbation with the same initial 6 and 6, unless it first collides with another mass element. 
The proof is trivial, following at once from equation (|17|) combined with the solution of the 
spherical model. 

Collision with another mass element inhibits collapse only if the matter is coUisional, 
in which case a shock wave is formed and the matter is compressed a finite amount. If the 
matter is collisionless, intersection of mass elements always increases 6 so that the collapse 
of each individual stream is speeded up. 

It is worth noting that the collapse theorem really makes two statements. First, any 
growing-mode perturbation with Si > (3/5) (fi^"^ — 1) collapses (or collides with another 
element) in a finite time. Second, shear speeds up the collapse. These two facts suggest 
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two corollaries. First, the collapse condition 5i > (3/5) — 1) may be too weak in the 
presence of shear: initially underdense mass elements may collapse. Second, because the 
nonlinear density growth rate depends on the shear, the initial shear affects the collapse 
time. Therefore, initial density maxima do not strictly identify the mass elements to collapse 
first. To determine how shear actually affects the collapse it is necessary to integrate all 
the evolution equations together, which we do in the next section under the assumpation of 
vanishing Hij. 

The collapse theorem applies regardless of the spatial geometry of the initial 
perturbation; it is based only on the trajectories of infinitesimally nearby particles. The 
theorem and its proof using the Raychaudhuri and mass conservation equations is similar to 
proposition 4.4.1 of Hawking & Ellis (1973), which underlies one of Hawking's singularity 
theorems. However, in the present case no spacetime singularity occurs (i.e., the spacetime 
remains geodesically complete) despite the collapse of matter to infinite density because 
we have assumed that the fluctuations are much smaller than the Hubble distance so that 
the Newtonian gravitational potential (eq. fl^), and therefore the metric perturbations, 
remains small. 



4. General Solutions for Hij = 

We have shown that the nonlinear evolution of density perturbations depends on the 
shear and tides. To describe the evolution fully we must integrate the equations of motion 
for these quantities in addition to the density and expansion scalar. We shall make the 
simplifying assumption that the magnetic part of the Weyl tensor, Hij, may be neglected. 



4.1. Overdense Perturbations 

Nonzero shear causes overdense perturbations to collapse sooner than the bound given 
by the collapse theorem. To find out how much sooner, we have integrated equations (^), 
(I), and ([TT|)-(p!^ numerically for Q = 1 with initial conditions given by equations (16). 
The expansion factor at collapse Oc was evaluated as a function of the initial tide parameters 
(eo, ao) for 6o = 1. Using equations (16) in the linear regime the results may be scaled to 
any 6o. Note that 6o refers to the linear amplitude at a = 1; the actual integrations began 
at a -C 1. 
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The results are shown in Figure 1. As expected, Oc decreases with increasing initial 
shear (or tide), verifying the collapse theorem. Noticing that the contours are approximately 
elliptical, we fit the following relation to our numerical results: 

a;i(5o>0,eo,ao) = Ci5o + C^2eo (1-0.2 cos ao) , = 1.25(1 - Ci) . (19a) 

The rms and maximum relative errors for in this fit are 1.5 and 3.6 percent, respectively, 
over the region shown in Figure 1. The fit fails to reproduce the asymmetry of the contours 
at large radii. An alternative fit, better for large eo, is 

a-\6^ > 0, eo, ao) =C,6o + eo (C3 + C, s + C, s"" + s^) (19b) 

where s = (1 - cos aoY^^ and C3 = 0.36950, C4 = -0.04219, C5 = 0.26710, and 
Cq = —0.07404. The rms and maximum relative errors are now 1.2 and 4 percent, 
respectively. 

Plane-parallel perturbations, with eo = So and cto = 0, represent a special case of 
collapse with a simple exact solution to our nonlinear equations: 6{t) = 5//(l — Si), 
9 = a = — 61), e = 61, and a = P = 0, where 5z(r) is the linear solution for the 

growing mode. This solution, valid for 5/ < 1, is identical to the result obtained by applying 
the standard Zel'dovich (1970) approximation ( [Matarrese et al. 1993| ). Croudace et al. 



(1993) have shown that the fully relativistic solution of this problem is given by the Szekeres 
(1975) metric. 

Zel'dovich (1970) suggested that, owing to kinematical effects, cold gravitational 
collapse would generically approach the plane-parallel solution. Surprisingly, we find this 
conclusion to be invalid for the solutions we have obtained. As 6 —>■ 00, the solutions 
generically approach a = 9 ^ —00 as in the Zel'dovich solution but with a + 2P = 3n and 
e diverging (roughly as 6) rather than remaining finite. Over much of the range of initial 
parameters we find a ~ tt at the moment of collapse, corresponding to a velocity gradient 
tensor with one positive and two negative (and equal) eigenvalues in the proportions 
1 : (—2 + a') : (—2 — a') where a' = (vr — a)/\/3 -C 1. In other words, the final state is 
two-dimensional collapse into a spindle or filament with expansion along the third axis, the 
opposite of what Zel'dovich predicted. This type of final configuration is achieved even if 
the behavior in the linear and weakly nonlinear regime is nearly one-dimensional, as we will 
see in §4.3 below. The final state is purely one-dimensional only if ao = 0. 

Croudace et al. (1993) noted the instability of the plane-parallel solution and suggested 
that it was due to the neglect of the magnetic part of the Weyl tensor. This is a possibility we 
cannot disprove because we have neglected Hij. However, we can determine why filaments 
form in this limit: the prolate configuration is favored by the nonlinear tide-shear coupling 
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term in equation (7), which is stabihzing for one- dimensional (pancake-hke) configurations 
with a ^ f3 ^ but destabihzing for two-dimensional (filamentary) configurations with 
a ^ P ^ 71. For a = P = 0, Eij has eigenvalues with signature (—,—,+) and aij has 
(+, +, — ) (assuming e > and o" < 0, the generic signs for collapse). The eigenvalues of 
their product (after subtracting the trace in eq. 0]) have signature (+, +, — ), so that this 
term, opposite in sign to Eij, retards the growth of tides. If a = /3 = vr, all of the signatures 
are reversed except for the tide-shear coupling term, which therefore becomes destabilizing. 
This behavior is seen most clearly in equation (|l^). It is possible that the magnetic term in 
equation (0) would, under some circumstances, counter the tide-shear coupling. However, 
this cannot be decided without further investigation. In any case, the results obtained here 
warn one not to take the generality of the Zel'dovich pancakes for granted given the neglect 
of nonlinear effects in the kinematical treatment. 

A complete physical explanation will have to await a Newtonian derivation of equation 
(0), on which this result hangs. However, the result is not surprising when one recalls 
that the gravitational field remains finite at a planar distribution of mass while it diverges 
at a linear distribution. While it is true that the divergence is greatest for a point-like 
(spherical) configuration, this state is never reached unless the shear vanishes. As we have 
seen, shear accelerates collapse, so that the spherical end state is disfavored. Collapsing 
mass elements find the next best solution, a prolate spindle. 

It is now easy to understand the asymmetry of the contours in Figure 1. Initially 
prolate configurations (with cosao < 0) collapse sooner than initially oblate configurations 
because the shear and tides are already oriented in the most unstable configuration. 

An important consequence of our result is that initial density maxima do not 
correspond, in general, to minima of the collapse time. Equations (19) show that the 
initial tide (or equivalently, shear) is of comparable importance to the initial density in 
determining the collapse time of a mass element. 



4.2. Underdense Perturbations 

Given the importance of shear and tides in accelerating the collapse of overdense 
perturbations, it is natural to inquire whether they can cause the collapse of underdense 
perturbations. To answer this question we integrated the equations of motion with Sq = —1. 
As Figure 2 shows, we found indeed that sufficiently large initial tide can turn an expanding 
underdense perturbation into a collapsing overdense one. Because the collapse time can 
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diverge (for sufficiently small initial tide), we graph a~^{So — —l,eo,a). This quantity 
vanishes at the innermost contour; configurations with smaller tide never collapse before 
trajectories intersect. Including the dependence on Sq, the numerical results are fit by 

3 / 1 \ 

a~^{6o < 0,eo,Q!o) = ^ C'i(eo + 25o) + ( C? - -^Cij (eocosao + 2So) + 

r 1 1/2 

Cy (eo cos Q!o + 2(5o)^ + 3eQ sin^ ao , (20) 

with C7 ~ 0.0676. The rms and maximum absolute errors in this fit for are 0.021 and 
0.046 over the range shown in Figure 2. The fit reproduces the kink at cq = —25o and 
tto = 0, but the location of the = contour does not match the numerical results 
exactly. Note that a negative implies that collapse does not occur. 

By comparing Figures 1 and 2 one may notice that for large cq the contours of constant 
Qc are very similar for do < and Sq > 0. In fact, excluding the kink region mentioned 
above, equation (19b), provides an excellent fit to both cases! Not only can negative density 
perturbations collapse in an Einstein-de Sitter universe, with large initial tide (eo ^ |5o|) 
they collapse at nearly the same rate as positive density perturbations. 

If the initial tide is small enough then a perturbation does not collapse before it 
intersects another mass element. In this case the mass drains out of the initial perturbation 
because its expansion rate always exceeds the Hubble expansion rate {9 > 0). At late times 
these solutions approach the expanding spherical void solution (Bertschinger 1985) with 
1 + S (X T"^, 9 — ^a/a, a (X — r^^lnr, and e = constant. The velocity field becomes radial 
{a/ 9 — > 0) despite the nonspherical gravitational potential {e/5 remains finite) because 
the net gravitational deceleration (proportional to 1 + 5) is so small that the perturbation 
becomes freely expanding, retaining thereafter whatever shape it had prior to entering the 
homologous expansion phase. 



4.3. General Perturbations 



It is worthwhile to unify the initial density and tide to obtain a representation of the 
solutions that is valid for any initial condition. We define a net perturbation A and an 
angle 7 so that the linear density perturbation and tide are 

5o = Ao COS70 , eo = ^ sin7o . (21) 

V5 
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The variables (Ao,7o, ao) are entirely equivalent to the eigenvalues of VjVj0 (eq. [|15|] ). 
The ranges of the new variables are Aq G [0, oo) and 70 G [0,7r]. The square root factor 
is introduced to diagonalize the covariance matrix of the variables (Ao,7o,ao) for a 
homogeneous and isotropic random process. 

If the linear density field So{x) is a gaussian random field, applying the results of 
Doroshkevich (1970) we find the joint probability distribution 

/ Aq \ Aq sin'^7o sinao 
dP(Ao, 7o, ao) = exp -— ^ ^ dAo d7o dao , (22) 

where as is the standard deviation of 6q. Equation (^2]) is reminiscent of a three-dimensional 
gaussian distribution in spherical coordinates aside from three extra powers of Aosin7o. 
These factors arise because the gravity gradient tensor is specified by three distinct diagonal 
entries of Eij (and angles giving the orientation of the eigenvectors which have been 
integrated out) in addition to 6. As a result, the distribution of Aq is shifted to larger 
values compared with a gaussian, while cos 70 is much more likely to be small than is cosao, 
which is uniformly distributed in [—1, 1]. We therefore change variables from 70 to 

16 no 2 4 1 

^0 = 1-7;- / sm 7 d7 = 1 7o + — sm 270 - — sm47o . (23) 

Sir Jo TT Sit dtt 

This variable is uniformly distributed (the measure dP depends on rjQ only through drjo) 

for a homogeneous and isotropic random process, ranging from rjo = ~1 (70 = tt) for 

tide-free negative density perturbations to r^o = +1 (7o = 0) for tide-free positive density 

perturbations. The distribution of the remaining variable, Aq, is a simple combination of 

gamma distributions of the argument of the exponential in equation ([2^). 

Figure 3 shows the inverse collapse factor (CiAo function of the two variables 

cosao and rjo, assuming irrotational growing-mode perturbations in an Einstein-de Sitter 
universe. Because these variables are independently uniformly distributed, this diagram 
may be viewed as an equal-area representation of the general case; all points in the plane 
are equally likely for homogeneous and isotropic random initial conditions. Using this 
interpretation and computing the fraction of the area below the = contour, we 
conclude that a randomly chosen mass element has a probability of 0.780 to collapse, 
implying that 56% of the underdense perturbations (and 100% of the overdense ones) will 
collapse. While the exact numerical values depend on the assumption that Hij = 0, the 
conclusion that some underdense regions must collapse is valid in general because of the 
effects of shear. Note that orbit-crossing only hastens collapse since it increases the density 
and therefore the gravitational focussing of trajectories. 

The figure presents our results in a rather different way than Figures 1 and 2. One's 
first impression may be that for a given Aq the fastest collapse occurs for the most spherical 
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initial perturbations. However, there are two flaws in this conclusion. First, although 
increasing rjo means increasing ratio of density to tide, this does not contradict our earlier 
conclusion that for fixed 5o, the slowest collapse occurs for spherical perturbations, because 
constant So and decreasing tide is not the same as constant Aq and increasing ratio of tide 
to density. One must be careful to specify what quantities are being held constant. The 
second problem is that we find a narrow region of large r/o where the (Aq ag)^^ exceeds the 
spherical collapse value Ci; Figure 3 shows the contour (CiAoac)"^ = 1.05. The contours 
double back at the top of this region so that (CiAq Oc)"^ = 1 on the fine r]o — 1. 

Although the collapse time appears to be a relatively simple function of cosao and tjq, 
we do not have yet an accurate fitting formula. Before investing much effort in this venture 
we should relax the assumption Hij = 0, a task we leave for future work. 

If we crudely approximate (CiAoflc)^^ ~ l^o + fj then local minima of the collapse 
time correspond to local maxima of Ao(|?7o + §)• It would be interesting to see in a 
realization of a gaussian random field how this quantity differs from ^o, whose maxima are 
widely (and erroneously) considered to be the first points to collapse. 

Figure 4 shows the final shape of the collapsed object, as represented by cos a, as a 
function of the initial parameters. Recalling that cos a = —1 for a two-dimensional (prolate) 
collapse while cos a — +1 for a one-dimensional (oblate) collapse, we see that, aside from 
regions with initially small shear (770 ~ ±1) and oblate shape (cosao ~ 1), most of the 
parameter space results in collapses that are nearly two-dimensional. In the shaded region 
the eigenvalues of the velocity gradient tensor deviate by less than 10 percent from the 
ratios 1 : —2 : —2 for a filament collapsing along two dimensions. This result confirms and 
extends what we found in §4.1 above. 



5. Conclusions 

In this paper wc have developed and applied a Lagrangian fluid approach to studying 
nonlinear gravitational instability. This approach, pioneered by Matarrese et al. (1993), 
has several advantages over the usual Eulerian description on one hand and the Zel'dovich 
approximation on the other hand. The flrst advantage is that we follow the fluid variables 
associated with a given mass element, which is what we want to do if we are to track the 
formation of objects by gravitational instability. Second, the Lagrangian fluid flow approach 
highlights the important physical role played by the shear and tides, quantities that are 
not so apparent in the Eulerian treatment and are generally neglected in the Zel'dovich 
approximation. 
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Going further requires making some assumption about the magnetic part of the 
Weyl tensor H^j, a quantity having no clear Newtonian interpretation but which may, 
nevertheless, be present in the Newtonian limit. If it is negligible, as we assumed in §4, then 
the Lagrangian approach has several additional advantages. The first and clearest is that 
the pressureless, irrotational flow problem reduces to a small set of six ordinary differential 
equations for, effectively, the eigenvalues of the velocity and gravity gradient tensors. These 
equations are easy to integrate to obtain accurate numerical solutions. Second, the method 
is exact (subject to the above-mentioned assumptions) until orbit-crossing. The Zel'dovich 
approximation is exact for one-dimensional perturbations. The Lagrangian method is exact 
for cylindrically and spherically symmetric perturbations also (in these cases Hij vanishes 
by symmetry). Third, this method is completely local in Lagrangian space (if Hij = 0) 
until trajectories intersect, so that one can follow the evolution of different mass elements 
independently up to orbit-crossing. The Poisson equation needs to be solved only once at 
the beginning to give the initial tidal field. 

These advantages come with a price: the Lagrangian fluid method does not keep track 
of particle trajectories. Thus, although one can compute how the density of a given mass 
element changes with time, one cannot know its position without separately integrating the 
equations of motion for the trajectory. The Lagrangian fluid equations could be integrated 
together with the trajectories in order to supplement standard N-body simulation methods. 
Matarrese et al. (1993) have demonstrated a more elegant method in which only differences 
in positions are integrated. With standard N-body techniques one can compute the density 
and shear only by averaging over volumes large enough to contain several particles. With 
the fluid flow approach the density and other quantities would be evolved for each particle, 
although the densities would have to be coarse-grained in regions of multi-streaming. 

In this paper we have stated and proven a theorem for pressureless, irrotational 
self-gravitating flows: For a given initial 5 and 5, a uniform spherical perturbation (i.e., one 
with zero shear) collapses (i.e., reaches infinite density) more slowly than any other shape. 
The proof requires no assumptions concerning the evolution of shear or tides or the magnetic 
part of the Weyl tensor. This theorem applies to the evolution of individual mass elements 
rather than to the global behavior of an extended mass distribution, for which spherical 
symmetry may still be preferred as a state of lowest energy. Moreover, if one defines 
coUapse by the requirement that all three axes collapse for an ellipsoidal perturbation, 
then the spherical tophat collapses most rapidly (White & Silk 1979). However, for many 
applications one wants to know when density singularities first form, in which case our 
collapse theorem is applicable. 

The collapse theorem contradicts standard lore. One's intuition suggests that uniform 
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spherical collapse ought to be most efficient because the gravitational field of a point mass 
is stronger than that for linear or planar distributions. However, this intuition neglects the 
very important role played by shear, which always acts to increase the rate of growth of 
fluid convergence —V ■ v following a given fluid element. The effect of shear is not apparent 
in the Euler equations until one writes the time derivatives following a given fluid element. 

Although exact analytical results are not currently available, we have numerically 
integrated the Lagrangian fluid equations assuming Hy = to determine the actual collapse 
time for perturbations with shear and to investigate the nature of the collapse. We found 
three interesting things. First, underdense regions can collapse if the initial shear is large 
enough. One might think this to be unlikely, but we showed that 56% of underdense 
perturbations (assuming a homogeneous and isotropic random process) will collapse if not 
prevented by shell-crossing. The exact percentage may differ if we relax the condition 
Hij = 0. 

Second, local minima of the collapse time do not correspond to local maxima of the 
initial density. Instead, they correspond to maxima of a combination of the eigenvalues of 
the initial gravity gradient (or velocity gradient) tensor. Although we have not obtained 
an analytical formula for this linear combination, we have given a crude approximation in 
terms of the variables Aq and r]o deflned in equations (^Tj) and (p3|). 

Third, we showed that when Hij = 0, the collapse is generically two-dimensional, 
in contrast with the prediction of Zel'dovich (1970) that gravitational collapse is one- 
dimensional. Zel'dovich did not consider the evolution of tides and he gave only a 
kinematical extrapolation of the linear evolution of trajectories. We flnd that, even if the 
initial collapse is nearly one-dimensional, nonlinear coupling between the tide and shear 
tends to drive the shape of the collapsing mass into a fllament rather than a pancake. We 
gave a heuristic explanation for this, noting that the gravitational fleld of a fllament is 
stronger than that of a pancake. Although the fleld of a point mass is still stronger, nonzero 
shear prevents uniform spherical collapse. 

Numerical simulations give conflicting answers to the question of whether cold 
gravitational collapse is one- or two-dimensional. Loeb & Rasio (1994) flnd that shear 
can lead to a fllamentary collapse for irrotational perturbations, but White (1993, private 
communication) found that a uniformly expanding triaxial ellipsoid collapses along one 
dimension flrst, in good agreement with the approximate analytic theory of White & Silk 
(1979) and in disagreement with the results obtained assuming Hij = 0. It appears that 
can cause a pancake-like collapse to occur in some circumstances. Giving Hij a Newtonian 
interpretation and determining the conditions under which it is signiflcant are important 
tasks for the future. 
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Our results have several important implications for the formation of galaxies and 
large-scale structure by gravitational instability. First, the widely-used spherical tophat 
collapse model underestimates the rate of gravitational collapse. A good deal of analytic 
theory in large scale structure has been based on the spherical collapse model, including the 
Press & Schechter (1974) theory for the mass function of collapsed objects. In a numerical 
simulation Katz et al. (1993b) noticed that objects with large initial shear collapse sooner as 
we would predict. Hoffman (1986) first made this point using the Zel'dovich approximation. 
Second, the two-dimensional character of at least some collapses means that galaxies may 
be expected to form along filaments. This point was also noted by Katz et al. Moreover, 
stronger initial tides for a given overdensity imply stronger prolateness of collapsing dark 
halos, leading to strong influences on the halo shape after virialization ( |Dubinski 1992| ). 
Third, initial density maxima are not the first objects to collapse. Indeed, Katz, Quinn, & 
Gelb (1993a) found from numerical simulations that many of the first objects to collapse 
in a cold dark matter model have no initial density maxima associated with them, calling 
into question the widely-used peak-bias model for galaxy formation. Fourth, matter in 
underdense regions also collapses unless the initial shear is too low. Thus, it is not clear that 
one expects galaxy formation to differ much in voids and in moderate-density environments. 

In future work we will investigate these implications. We will also test the results using 
N-body simulations and, conversely, test the N-body method by its ability to reproduce 
our theoretical results. A final important area of study concerns the evolution of the 
Weyl tensor. Although we have borrowed results from general relativity, the equations of 
motion for Eij and Hij should be derivable from Newton's laws. Conversely, the Newtonian 
description should help to clarify the proper relativistic treatment after orbits cross as well 
as the interpretation of Hij. 

We wish to thank Andrew Hamilton, Rien van de Weygaert, Simon White, and Fred 
Rasio for useful discussions. This work was supported by NSF grant AST90-01762. 
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6. Figure Captions 

Fig. 1: Contours of constant collapse time, as indicated by the cosmic expansion factor 
Oc, versus initial tidal field parameters for linear density perturbation Sq — +1. The 
light (heavy) contours are spaced by 0.1 (0.5), with the central contour Oc = 1.6 and the 
outermost contour Oc = 0.4. The maximum collapse time, Oc = 1.686, occurs at the center 
corresponding to zero shear. 

Fig. 2: Contours of constant inverse collapse factor for initial negative density 
perturbations with 6o — —1. The light (heavy) contours are spaced by 0.1 (0.5), with the 
innermost contour — and the outermost one — 1.8. Initial perturbations in the 
central region do not collapse, while perturbations with > Ci — 0.593 (oc < 1.686) 
collapse faster than spherical positive density perturbations of the same \Sq\. 

Fig. 3: Contours of constant inverse collapse factor a~^, normalized to the spherical value 
Ci, for perturbations with Aq = 1 (see eq. [21] for the definition). The shape of the initial 
tidal tensor is specified by cos«0) with cos«o = ~1 (+1) for prolate (oblate) perturbations. 
The parameter rjo defines the ratio of density to tide, with rjo = +1 (—1) for spherical 
positive (negative) density perturbations and i]q = for pure tidal perturbations with 
6o = 0. The light (heavy) contours are spaced by 0.1 (0.5), aside from the uppermost 
contour, which is — 1.05 Ci. The lowest contour is — 0, below which perturbations 
do not collapse. 

Fig. 4: Contours of constant final shear shape cos a versus the same initial tidal parameters 
shown in Figure 3. Starting in the upper right corner and moving to the lower left, the 
contours are (heavy), —.1, —.2, . . ., —.9 (heavy), —.95, —.98, —.99, —.999 (heavy), and 
thereafter they reverse this sequence. Over most of this plane the final collapse is strongly 
prolate (cos a —1). The shaded region shows cosa < —0.95. 



